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Abstract. The phase space trajectories of many body systems charater- 
istic of simple fluids are highly unstable. We quantify this instability by a 
set of Lyapunov exponents, which are the rates of exponential divergence, 
or convergence, of initial (infinitesimal) perturbations along carefully se- 
lected directions in phase space. It is demonstrated that the perturbation 
associated with the maximum Lyapunov exponent is localized in space. 
This localization persists in the large-particle limit, regardless of the in- 
teraction potential. The perturbations belonging to the smallest positive 
exponents, however, are sensitive to the potential. For hard particles they 
form well-defined long-wavelength modes. The modes could not be ob- 
served for systems interacting with a soft potential due to surprisingly 
large fluctuations of the local time-dependent exponents. 



1 Lyapunov spectra 

Recently, molecular dynamics simulations have been used to study many body 
systems representing simple fluids or solids from the point of view of dynamical 
systems theory. Due to the convex dispersive surface of the atoms, the phase 
trajectory of such systems is highly unstable and leads to an exponential growth, 
or decay, of small (infinitesimal) perturbations of an initial state along specified 
directions in phase space. This so-called Lyapunov instability is described by a 
set of rate constants, the Lyapunov exponents {Xi,l — 1, . . . , D}, to which we 
refer as the Lyapunov spectrum. Conventionally, the exponents are taken to be 
ordered by size, A; > A; + i. There are altogether D = 2dN exponents, where d 
and D denote the dimensions of space and of phase space, respectively, and N 
is the number of particles. For fluids in nonequilibrium steady states close links 
between the Lyapunov spectrum and macroscopic dynamical properties, such 
as transport coefficients, irreversible entropy production, and the Second Law 
of thermodynamics, have been established ■ This important result 

provided the motivation for us to examine the spatial structure of the perturbed 
states associated with the various exponents. Here we present some of our results 
for two simple many-body systems representing dense two-dimensional fluids 
in thermodynamic equilibrium. The first model consists of N hard disks (HD) 
interacting with hard elastic collisions, the second of N soft disks interacting 
with a purely repulsive Weeks-Chandler-Anderson (WCA) potential. 
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The instantaneous state of a planar particle system is given by the AN- 
dimensional phase space vector T — {rj,pj,;i = 1,...,N}, where and p; 
denote the respective position and linear momentum of molecule i. An infinites- 
imal perturbation 5T = {Sri,5pi;i — l,...,N} evolves according to motion 
equations obtained by linearizing the equations of motion for T(t). For ergodic 
systems there exist D = AN orthonormal initial vectors { (51^(0); 1 = 1,..., 4N} 
in tangent space, such that the Lyapunov exponents 



exist and are independent of the initial state [|7fl8|] . Geometrically, the Lyapunov 
spectrum describes the stretching and contraction along linearly-independent 
phase space directions of an infinitesimal hypersphere co-moving with the flow. 
For equilibrium systems the symplectic nature of the motion equations assures 
the conjugate pairing rule to hold |J: the exponents appear in pairs with a 
vanishing pair sum, A;+A4tv +1 _; = 0. Thus, only the positive half of the spectrum 
{Ak;<2At} needs to be calculated. The sum of all Lyapunov exponents vanishes, 
which according to Liouville's theorem is a manifestation of the conservation of 
phase volume by Hamiltonian systems. Six of the exponents, {A2jv-2<z<2JV+3}, 
always vanish as a consequence of the conservation of energy, momentum, and 
center of mass, and of the non-exponential time evolution of a perturbation 
vector parallel to the phase flow. 

For the computation of a complete spectrum a variant of the classical algo- 
rithm by Benettin et al. and Shimada et al. is commonly used j^,[nj • It requires 
following the time evolution of the reference trajectory and of an orthonormal 
set of tangent vectors {8Ti{t);l — 1, . . . ,4iV}, where the latter is periodically 
re-orthonormalized with a Gram-Schmidt (GS) procedure after consecutive time 
intervals Atas- The Lyapunov exponents are determined from the time- averaged 
renormalization factors. For the hard disk systems the free evolution of the par- 
ticles is interrupted by hard elastic collisions and a linearized collision map needs 
to be calculated. An algorithm based on this idea has been developed by Del- 
lago et al. |Lq ], Although we make use of the conjugate pairing symmetry and 
compute only the positive branch of the spectrum, we are presently restricted 
to about 1000 particles by our available computer resources. 

For our numerical work reduced units are used. In the case of the Weeks- 
Chandler-Anderson interaction potential, 



r 
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the particle mass m, the particle diameter a, and the time (ma 2 / e) 1 / 2 are unity. 
In this paper we restrict our discussion to a thermodynamic state with a total 
energy per particle, E/N, also equal unity. For the hard-disk fluid, (Nmu 2 / K) 1 / 2 
is taken as the unit of time, where K is the total kinetic energy, equal to the 
total energy E of the system. There is no potential energy in this case. For the 
reduced temperature we have T = K/N, where Boltzmann's constant is also 
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unity. In the following, Lyapunov exponents for the two model fluids will be 
compared for equal temperatures (and not for equal total energy) . This requires 
a rescaling of the hard-disk exponents by a factor of \J {K/N)wca/(K/N)hd to 
account for the difference in temperature. All our simulations are for a reduced 
density p = N/V — 0.7, where the simulation box is a square with a volume 
V = L 2 and a side length L. Periodic boundaries are used throughout. 




I/2N 

Fig. 1. Lyapunov spectrum of a dense two-dimensional fluid consisting of N = 
100 particles at a density p = 0.7 and a temperature T = 0.75. The label 
WCA refers to a smooth Weeks-Chandler- Anderson interaction potential, and 
HD indicates the hard disk system. 



As an example, in Fig. |l| we compare the Lyapunov spectrum of a 100- 
particle WCA fluid to a spectrum for an analogous hard disk system at the 
same temperature (T — 0.75) and density (p = 0.7). A reduced index 1/2N 
is used on the abscissa. It is surprising that the two spectra differ so much in 
shape and size. The difference persists in the thermodynamic limit. The step-like 
structure displayed by the hard disk spectrum for 1/2N close to 1 is an indication 
of a coherent wave-like shape of the associated perturbations, to which we refer 
as Lyapunov modes fl3|,|l4| . We defer the discussion of these modes to Section 
4. 



4 



2 Time-dependent local exponents 

We infer from Equ. (|l|) that the Lyapunov exponents are time averages over an 
(infinitely) long trajectory and, therefore, are global properties of the system. 
They can be rewritten as 



r — >oo 



A; = lim / A',(T(t))<tt = (A'i), (3) 



where the (implicitly) time-dependent function A';(r) depends on the state r(i) 
which the system occupies in phase space at time t. Thus, A'z(r) is called a local 
Lyapunov exponent. It may be estimated from 

,, fr ,,u i , mjTjt + At GS )\ 

where t and t+Atcs refer to times immediately after consecutive Gram-Schmidt 
re-orthonormalization steps. Its time average along a trajectory is denoted by 
(• • •) and gives the global exponent A;. 
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Fig. 2. Second moment of all local Lyapunov exponents for a planar 16-particle 
system. The labels WCA and HD refer to the Weeks-Chandler-Anderson and 
hard disk models, respectively. The re-orthonormalization time Abas = 0.075, 
the temperature T = 0.75, and the density p = 0.7. 
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The local exponents fluctuate considerably along a trajectory. This is demon- 
strated in Fig. [| where we have plotted the second momens (A' 2 ) as a function 
of the Lyapunov index 1,1 < I < 47V , for a system of 16 particles, both for the 
WCA and HD models. I = 1 refers to the maximum, and I — 64 to the most- 
negative exponent. The points for 30 < I < 35 correspond to the 6 vanishing 
exponents and are not shown. We infer from this figure that for the soft WCA 
particles the fluctuations of the local exponents become very large for the Lya- 
punov exponents describing relatively-weak instabilities with near-zero growth 
rates, / — > 2N. For the hard disk system, however, the relative importance of the 
fluctuations becomes minimal for the same exponents. We shall return to this 
point in Section 4. 

We note that the computation of the second moments (A' 2 ) for the hard 
disk system requires some care. Due to the hard core collisions they depend 
strongly on Atcs- The variance of the fluctuating local exponents, (A' 2 } — (A';) 2 , 
varies with l/Atcs f° r small Atcs, a s is demonstrated in Fig. @ for I = 1. 
However, the shape of the fluctuation spectrum is hardly affected by the size of 
the renormalization interval. 
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Fig. 3. Mean-square fluctuations of the local exponent A'i, (A' 2 } — (A'i) 2 , times 
Atos, as a function of the re-orthonormalization interval Atcs- The system 
consists of 64 hard disks. The temperature T — K/64 — 1, and the density 
p = 0.7. 
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3 The maximum exponent 

The maximum Lyapunov exponent is the rate constant for the fastest growth 
of a phase space perturbation in a system. Thus, it is dominated by the fastest 
dynamical events in the fluid, binary collisions. There is strong numerical evi- 
dence for the existence of the thermodynamic limit {A^, V — > oo, TV/V^constant} 
for Ai and, hence, for the whole spectrum JT2|Jl4| . Furthermore, the pertur- 
bation belonging to Ai is strongly localized in space Jl6| , |l7 14 . This may be 



demonstrated by projecting the tangent vector STi onto the four-dimensional 
subspaces spanned by the perturbation components of the individual particles. 
The squared norm of this projection, Sf(t) = (8ri)f + (Spi)f , indicates how active 
a particle i is engaged in the growth process of the pertubation associated with 

Ai pm 




Fig. 4. The surface depicts the instantaneous contributions of the individual 
particles, plotted at the position of the particles in the simulation cell, to the 
squared norm of the perturbation vector associated with the maximum exponent 
Ai. The system consists of 100 hard disks. The temperature T = K/N = 1, and 
the density p = 0.7. 



In Fig. U 5f(t) is plotted along the vertical axis for all particles of a hard 
disk system at the respective positions (xi,yi) of the disks in space, and the 
ensuing surface is interpolated over a periodic grid covering the simulation box. 
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The figure refers to an instantaneous condition for a well-relaxed system, and no 
averaging is involved. A strong localization of the active particles is observed at 
any instant of time. Similar, albeit slightly broader peaks are observed for the 
WCA system or other soft disk potentials. 

This localization is a consequence of two mechanisms: firstly, after a colli- 
sion the delta- vector components of two colliding molecules are linear functions 
of their pre-collision values and have only a chance of further growth if their 
values before the collision were already far above average. Secondly, each renor- 
malization step tends to reduce the (already small) components of the other 
non-colliding particles even further. Thus, the competition for maximum growth 
of tangent vector components favors the collision pair with the largest compo- 
nents. The active zone moves in space in a diffusive manner. 

The localization also persists in the thermodynamic limit. To show this we fol- 
low Milanovic et al. (m) and order all squared components [<5ri]|; j = 1, . . . ,4N 
of the perturbation vector 8Ti according to size. By adding them up, starting 
with the largest, we determine the smallest number A of terms, required for the 
sum to exceed a threshold 0. Then, C\ : e = A/ AN may be taken as a relative 
measure for the number of components actively contributing to Ai; 

/4JVCi, e \ 

E l ST tfs)> [S^ > [6^ for i<j. (5) 

Obviously, Ci t \ = 1. 

In Fig. U C%,e is shown for — 0.90 as a function of the particle number 
N, both for the WCA fluid and for the hard disk system. It converges to zero if 
our data are extrapolated to the thermodynamic limit, N — > oo. This supports 
our assertion that in an infinite system only a vanishing part of the tangent- 
vector components (and, hence, of the particles) contributes significantly to the 
maximum Lyapunov exponent at any instant of time. 



4 Lyapunov modes 

We have already mentioned the appearance of a step-like structure in the Lya- 
punov sepctrum of the hard disk system for the positive exponents closest to 
zero, which was first observed in Ref. | fl2| . They are a consequence of coher- 
ent wave-like spatial patterns generated by the perturbation vector components 
associated with the individual particles [ ^3[|l^| . 

In Fig. ^this is visualized by plotting the perturbations in the x (bottom surface) 
and y directions (top surface), {Sxi, i — 1, . . . , N} and {Syt, i = 1, . . . , TV}, re- 
spectively, along the vertical axis at the instantaneous particle positions (xj, j/j) 
of all particles i. This figure depicts a transverse Lyapunov mode of type Ti, 
for which the perturbation components are perpendicular to the respective wave 
vectors with one wave length equal to L, the linear extension of the simulation 
box. The system consists of N = 1024 hard disks, and the perturbation vec- 
tor considered in the figure belongs to the smallest positive exponent A2045- An 
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Fig. 5. Dependence of the localization measure, Ci t @, on the number of particles 
N for the parturbation associated with the maximum Lyapunov exponent Ai. 
The threshold = 0.90. The soft and hard disk results are labeled by WCA and 
HD, respectively. 
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1=2045 




Fig. 6. Transverse Lyapunov mode of type T\ associated with the smallest pos- 
itive exponent A2045 for 1024 hard disks at a density p = 0.7. The temperature 
K/N = 1. 

Bottom: Plot of the tangent vector components Sx of the individual particles 
(along the vertical axis) for all particles at their instantaneous postions (x, y) in 
the simulation box (horizontal plane). Top: Analogous plot for Sy. 
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analogous plot of Sp x and 5p y for the same perturbation vector is identical in 
shape to that of dx and 5y in Fig. 6, with the same phases of the waves. This is a 
consequence of the fact that the perturbations are solutions of linear first-order 
equations instead of second order equations such as the wave equation. Further- 
more, the exponents for 2042 < I < 2045 are equal. This four-fold degeneracy of 
non-propagating transversal modes, and an analogous eight-fold degeneracy of 
propagating longitudinal modes, are responsible for a complicated step structure 
for I close to 2N, which has been studied in detail in Refs. Jl3| , ^9[ . 

The wave length of the modes and the value of the corresponding exponents 
are determined by the size L of the square simulation box. There is a kind of 
linear dispersion relation [[l3|Jl4| indicating that the smallest positive exponent 
is proportional to 1/L. This assures that for a square simulation box with aspect 
ratio 1 there is no positive lower bound for the positive exponents of a hard disk 
system in the thermodynamic limit [^4|Jl9[| . 

So far, our discussion of modes applies only to hard disk fluids. In spite of a 
considerable computational effort we have not yet been able to indentify modes 
for two-dimensional fluid systems with a soft interaction potential such as WCA 
or similar potentials |po|] . This surprising result seems to be due to the very 
strong fluctuations of the local exponents discussed in Section 2. They tend to 
obscure any mode in the system in spite of considerable averaging and make a 
positive identification very difficult. 

We are grateful to Christoph Dellago, Robin Hirschl, Bill Hoover, and Ljubo- 
mir Milanovic for many illuminating discussions. This work was supported by the 
Austrian Fonds zur Forderung der wissenschaftlichen Forschung, Grants P 11428- 
PHY and P15348-PHY. 
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